Impurity corrections to the thermodynamics in spin chains using a transfer-matrix 

DMRG method 
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We use the density matrix renormalization group (DMRG) 
for transfer matrices to numerically calculate impurity correc- 
tions to thermodynamic properties. The method is applied 
to two impurity models in the spin- 1/2 chain, namely a weak 
link in the chain and an external impurity spin. The numer- 
ical analysis confirms the field theory calculations and gives 
new results for the crossover behavior. 



I. INTRODUCTION 

The study of quantum impurities remains a large part 
of condensed matter physics. The Kondo effect is maybe 
one of the most famous examples for impurity effects, but 
more recently much effort has been devoted to impurities 
in low-dimensional magnetic systemadn connection with 
high temperature superconductivityfcl. For the particu- 
lar case of impurities in quasi one-dimensional systems, 
much progress has been made with field theory descrip- 
tions e.g. for the Kondo model 2 !, quantum wiresa, and 
spin chainso. In those cases, the impurity behaves effec- 
tively as a boundary condition at low temperatures and 
the behavior can be described in terms of a renormal- 
ization crossover between fixed points as a function of 
temperatureB. 

Numerically this renormalization picture has been 
confirmedeTEl, but so far it was only possible to examine 
a limited number of energy eigenvalues in the spectrum. 
While some efforts have been made to extract thermo- 
dynamic properties from the energy spectrum directlyQ, 
such an approach is tedious and remains limited by finite 
system sizes. Monte Carlo simulations appear to be well 
suited for determining thermodynamic properties, but for 
the particular case of impurity properties it turns out to 
be difficult to accurately determine a correction which 
is of order 1/N, where N is the system size. We now 
apply the transfer matrix DMRG to impurity systems. 
This overcomes those problems by explicitly taking the 
thermodynamic limit N — ► oo, while still being able to 
probe impurity corrections and local properties at finite 



temperatures even for frustrated systems (which are not 
suitable for Monte Carlo simulations due to the minus 
sign problem). 

There are two separate impurity effects that we wish 
to address. The first is the impurity correction .Pimp to 
the total free energy of a one dimensional system 



Ft — 

± imp — 



lim (Ftotal 



NF pnie ) 



(1) 



where impure is the free energy per site for an infinite sys- 
tem without impurities. In other words, the impurity 
contribution is that part of the total free energy that 
does not scale with the system size N 
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Therefore, the impurity free energy Fi mp is directly pro- 
portional to the impurity density of the system, and it 
immediately determines the corresponding impurity spe- 
cific heat and impurity susceptibility, i.e. quantities that 
can be measured by experiments as a function of temper- 
ature and impurity density. Despite the obvious impor- 
tance of this impurity contribution we are not aware of 
any numerical studies that considered this quantity for 
any non-integrable impurity system. Traditional meth- 
ods would require an extensive finite size scaling analysis 
to track down the 1/N correction to the total free en- 
ergy per site, but our approach allows us to calculate 
-Pimp directly in the thermodynamic limit. We would 
like to point out that in other studies the response of an 
impurity spin to a local magnetic field is often termed 
"impurity susceptibility" , but we prefer to reserve this 
expression for the impurity contribution 
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where B is a global magnetic field on the total system. 

The second aspect of impurity effects are local proper- 
ties of individual sites near the impurity location, e.g. cor- 
relation functions and the response to a local magnetic 
field. Local properties have been the central-, uart of a 
number of works for many impurity modelsErt2l. Our 
approach is now able to calculate these impurity effects 
directly in the thermodynamic limit and we get quick 
and accurate results to extremely low temperatures. It 
turns out that the local impurity effects can be deter- 
mined much more accurately and to lower temperatures 
than the impurity contribution -Fi mp , which remains lim- 
ited by accuracy problems even with this method. 
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The Density _ . Matrix Renormal- 

ization Group (DMRG)EJ has had a tremendous success 
in describing low energy static properties of many one- 
dimensional (ID) quantum systems such as spin chains 
and electron systems. More recently many useful exten- 
sions to the DMRG have been developed. Nishino showed 
how to successfully apply the density matrix idea to two- 
dimensional (2D) classical systems by— determining the 
largest eigenvalue of a transfer matrixlij. Bursill, Xiang 
and Gehring have then shown that the same idea can be 
used to calculate thermodynamic properties of the quan- 
tum spin- 1/2 XY"-chaint2l. The method has later-been 
improved by Wang and XiangJiJ as well as ShibataEij and 
been applied to the anisotropic spin-1/2 Heisenberg chain 
with great success. In this paper we apply a generaliza- 
tion of the method to impurity systems, which is pre- 
sented in section [n| In Section III we study two different 
impurity models in the spin-1/2 chain and are able to 
confirm predictions from field theory calculations. Sec- 
tion IV concludes this work with a discussion of the re- 
sults and a critical analysis of accuracy and applicability 
to other systems. 



II. THE METHOD 

The method of the transfer matrix DMRG can in prin- 
ciple be applied to any one-dimensional system for which 
a transfer matrix can be defined. As a concrete exam- 
ple we will consider the antiferromagnetic spin-1/2 chain, 
since this model is well understood in terms of field theory 
treatments and has direct experimental relevance. The 
Heisenberg Hamiltonian can be written as 



A' 



H = hi 



Bi Sf 



(4) 



where Ji is the exchange coupling between sites i and i + 
1, and Bi is an external magnetic field in the z-direction 
at site i. Periodic boundary conditions, Sn + i = Si, are 
assumed. The partition function is defined by 



= h- C -0H _ 



Z = tr e 
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(5) 



where (3 — -^Kp and where we in the last step have par- 
titioned H into odd and even site terms, 



N/2 N/2 



(6) 



i=i 



A. The transfer matrix method 



The quantum transfer matrix for this system-is defined 
as usual via the Trotter-Suzuki decompositionlij 
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where M is the Trotter number. This expression ap- 
proximates the partition function up to an error of order 
(j3/M) 2 and becomes exact in the limit M — ► oo. By 
inserting a complete set of states between each of the ex- 
ponentials in Eq. (Q) and rearranging the resulting matrix 
elements, the partition function canJae written as a trace 
over a product of transfer matrices,E3 
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where Tm{2i — 1) is the 2 2M x 2 2M dimensional quantum 
transfer matrix from lattice site 2i — 1 to site 2i + 1. Note 
that Tm is in general not symmetric. However, if the two- 
site Hamiltonian h 2 i-i and h 2 i of Eq. (^) is invariant un- 
der the exchange 2i—l <-> 2i, and 2i <-» 2z+l respectively, 
as is the case unless we have applied a non-uniform mag- 
netic field, Tm is a product of two symmetric transfer 
matrices, one from site 2i — 1 to site 2i and the other 
from site 2i to 2i + 1. For a uniform system, the transfer 
matrix is independent of lattice site, Tm(2i — 1) = Tm, 
and the partition function is then given by 
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In the thermodynamic limit of the uniform system, Zm 
is given by 



N- 



lim Z M = X N/2 
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where A is the largest eigenvalue of Tm ■ 

The largest eigenvalue A can be found exactly only 
for small Trotter numbers M. As M increases, the di- 
mension of Tm grows exponentially, and we have to use 
some approximation technique to find A. Analogous to 
the case where the DMRG can be used to find a cer- 
tain eigenstate of a Hamiltonian as the number of lattice 
sites increase, we can use the DMRG to find the largest 
eigenvalue of Tm as the Trotter number M increase. The 
strategy is thus to start with a system block, T M ^ 21 and 
an environment block, T M ^ 2 , with a small M. The su- 
perblock transfer matrix, Tm, with Trotter number M, is 
constructed by "gluing" together the system block with 
the environment block. Periodic boundary conditions in 
the Trotter direction must be used. The reduced density 
matrix is constructed from the target state, i.e. the eigen- 
state of Tm with largest eigenvalue. Since the transfer 
matrix is non-symmetric, the left and right eigenvectors 
will not be complex conjugates of each other. 

A reduced density matrix for the system as part of the 
superblock can be constructed by taking a partiaL trace 
of T^J 2 over the environment degrees of freedomt-3, 



o - — tr T N/2 
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In the thermodynamic limit only the state with the 
largest eigenvalue will contribute, 



P 



N- 



(12) 



where \ip R ) and (ip L \ are the right and left eigenvectors 
of the superblock transfer matrix, Tm, corresponding to 
the largest eigenvalue, A. The matrix elements are given 
by Pi',i — J2j ^Vj^i,] > wnere % and j label the degrees of 
freedom of the system and the environment respectively 
and the target states are given by = J2i j ^j(*IO'l 

and|V^)=E lJ ^N)b'>- 

The left and right eigenvectors of the density matrix 
with the largest eigenvalues are then used to define the 
projection operators onto the truncated basis. After the 
first iterations we will keep m states for the system and 
the environment and the superblock will be Am 2 dimen- 
sional. More details on the transfer matrix DMRG al- 
gorithm for quantum systems have been presented in 
Refs. H-H 

Since p is non-symmetric it is not obvious that the 
eigenvalues are real, but because p represents a density 
matrix we expect it to be positive definite. However, be- 
cause of numerical inaccuracies complex eigenvalues tend 
to appear in the course of iterations. This is usually 
connected to level crossings in the eigenvalue spectra of 
the density matrix as M is increased. In addition, we 
also observe that just before the complex eigenvalues ap- 
pear, multiplet symmetries in the eigenvalue spectrum 
of p may split up (often due to "level repulsion" from 
lower lying states which were previously neglected). To 
overcome this problem it is important to keep all original 
symmetries (i.e. eigenvalue multiplets of p) even as M in- 
creases. For that purpose it is often necessary to increase 
the number of states m just before a multiplet tends to 
split up, which avoids the numerical error that leads to 
the symmetry breaking. In case of a level crossing of two 
multiplets that do not split up, complex eigenvalues may 
still appear and it is then possible to numerically trans- 
form the complex eigenstate pair into a real pair spanning 
the same space. Since the transformed real pair still is 
orthonormal to every other eigenvector, this transforma- 
tion does not cause any troubles for later iterations. With 
this method the eigenvalues stay complex until the two 
levels have crossed and moved enough apart, at which 
point the eigenvalues become real again. The renormal- 
ization procedure is therefore not as straightforward as 
in ordinary DMRG runs, since it is essential to track all 
eigenvalues of p and to dynamically adjust the number of 
states m, thereby sometimes repeating previous iteration 
steps. 



B. Impurities and local properties 

The renormalization scheme also allows the calculation 
of thermal expectation values of local operators. The 
magnetization of the spin at site i is determined by 



(S!) = - tr Sie-P*. 



(13) 



The operator Sf thus only has to be incorporated in the 
corresponding transfer matrix at site i. For the pure sys- 
tem, we then arrive at the formulaGJ 



m 1 



A 



(14) 



where T^(i) is defined similar to Tm(i) but with the 
additional operator Sf included in addition to the Boltz- 
mann weights. To measure the local bond energy, h(i) = 
Si ■ Sj.fi, a similar construction can be used. 

Let us now assume that the system has a single impu- 
rity. The systems we will study is the periodic spin-1/2 
chain with one weakened link or an external spin-1/2 
coupled to one spin in the chain 



Hi — Hg — SJ Sjv • Si 
H 2 = H + J' Sx ■ Sf, 



where 



JV-l 



Hq — ^ JSi ■ Si + i + JSn ■ Si 



(15) 
(16) 



(17) 



is the periodic chain and S/ is the spin operator of an 
external spin-1/2. The models are depicted in Fig. [l]. 
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FIG. 1. The impurity configurations considered in 
Eqs. ( |l5| ) and (|l6|). Note, that the left and right ends are 
joined with periodic boundary conditions. 

For systems with such a local impurity, which is con- 
tained within two neighboring links, only one of the 
Tjvf(i) in Eq. (|J) will differ from a common Tm 
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tr (t*/ 2 - 1 ^ 



- M 
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where Ti mp is the transfer matrix of the two links contain- 
ing the impurity and Tm is the transfer matrix describing 
the bulk. In the thermodynamic limit the partition func- 
tion will still be dominated by the largest eigenvalue of 
the "pure" transfer matrix. From Eq. (|l8| ) we have 



N 



lim Z M = \ N / 2 ~ 1 (4> L \T imp \4> R ), 



(19) 
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where A, (ip L | and \ip R ) all correspond to the pure system. 
The generalization of Eqs. (|l^ ) and (|l^) to impurity con- 
figurations ranging over more than two links is straight- 
forward. In this case more than one impurity transfer 
matrix has to be introduced and this could be used to 
study e.g. multiple impurities and impurity-impurity in- 
teractions. 

Let us define A; mp = (^ L \Ti mp \ijj R ). The total free 
energy of the system is then given by 



F 



1 



InZ 



lln (A^ 2 -^) 



N 



■In A 



1 ln Aimp 



2(i (3 A 

By comparing with Eq. (^|) we can retrieve the pure and 
impure parts 



F pure = -ilnA, F imp --iln^. 



(21) 



The impurity susceptibility can be found from the change 
of -Pimp in a small magnetic field from Eq. (|^) 



Ximp 



dB 



2 x imp ■ 
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Local properties such as the magnetization of the impu- 
rity spin can be determined by 



(^imp) 



(^|T- p |^) 



A, 



(23) 



The magnetization of spins close to the impurity is read- 
ily obtained by 



(S z ) 



^ L \Tll{T M ) x T lmv \^ R ) 

A AmiD 



(24) 



where 2x is the number of sites between the impurity 
and the spin of interest. Note that since a transfer ma- 
trix involves a total of three lattice sites, Tfjf can be 
constructed to measure the spin at any of these sites (or 
the mean value) . The actual site of the measurement in 
Eq. ( gjj ) is thus determined both by x and how T|| is 
set up. The expectation value in Eq. ( pi]) is most easily 
calculated by first computing the vectors (ip L \T^ and 
Tim P \ip R }, then acting with Tm on one of these states, 
and finally calculating the inner product of the resulting 
states. Eq. ( p4[ ) can be generalized to measure any equal 
time correlation function with or without an impurity, 
e.g. by replacing T imp by T|f. 

The reduced density matrix for the impurity system 
can be constructed by taking the thermodynamic limit 
of the impurity version of Eq. (ill), 



1 

Zm 



P = 7^ tr onv I T, 



N/2- 
M 



imp 



(25) 



with Zm as in Eq. (|l8|). In our calculations we have 
used the same density matrix as for the pure case, i.e. 
Eq. ([l2|) ; and we have found it to give good results. This 
form can most easily be motivated by writing Eq. (|2E 



on the form p = tr cnv (t^/ 4 L T imp T M / ' 1 ) /Z M . From 

a computational point of view, this method is also very 
convenient; by storing all target states and projection 
operators from the DMRG run for the pure system, all 
local impurities can be studied by simply using the same 
projection operators and target states. This makes sub- 
sequent DMRG runs for different impurity parameters 
very fast. 

There are also other choices of density matrices that 
can be made. The thermodynamic limit of Eq. ( p5|) can 
(20) a l so DC interpreted as p — * IV^KV^I^lmp/Aimp in which 



case the impurity transfer matrix would be taken into 
account in the density matrix. This would in some sense 
be analogous to including an operator different from 
the Hamiltonian in the density matrix of the ordinary 
"zero-temperature" DMRG, which is usually not neces- 
sary to measure e.g. correlation functions in the ordinary 
DMRG. This approach would destroy the computational 
advantage of using the pure density matrix, because the 
pure projection operators and target states could not be 
used but instead complete DMRG runs would have to be 
done for each impurity configuration and coupling. 



III. RESULTS 
A. Field theory predictions 

To make a meaningful analysis of the numerical re- 
sults, we first need to understand the spin-1/2 chain in 
the framework of the quantum field theory treatment. 
This turns out to give a good description of the impurity 
behavior in terms of a renormalization flow between fixed 
points, and we will be able to set up concrete expecta- 
tions for the impurity susceptibility in Eq. (||) as well as 
local properties. 

The effective low energy spectrum of the spin-1/2 chain 
is well described by a free boson Hamiltonian density 



n = ~ [(n,) 2 



(26) 



plus a marginal irrelevant operator cos v87T0 and other 
higher order operators which we have neglected. Here 
is the momentum variable conjugate to </>. In the long 
wave-length limit the spin operators can be expressed in 
terms of the bosonic fields using the notation of Ref. 



s; 



l 



2?r dx 



+ (— l)- J const. cosv27T0 



iV 27T0 



[const. cosV2~K<p + ( — l) J const.], (27) 



At t his p oint we can introduce the impurities in Eq. ( |l5| ) 
and <M) in a straightforward way as perturbations. The 
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field theoretical expressions for these perturbations can 
then be analyzed in terms of their leading scaling dimen- 
sions. Local perturbations with a scaling dimension of 
d > 1 are considered irrelevant, while perturbations with 
a small scaling dimension d < 1 are relevant and drive 
the system to a different fixed point. Hence, we can pre- 
dict a systematic renormalization flow towards or away 
the corresponding fixed point, respectively. 

Such an analysis has been made in Ref. ^ and the renor- 
malization flows have been confirmed by determining; the 
finite size corrections to the low energy spectrumo. In 
particular, a small weakening of one link in the chain 

Hi =H -SJS N -Si (28) 

has been found to be a relevant perturbation described by 
the operator sin-\/27r</> with scaling dimension d = 1/2, 
so that the periodic chain (6J = 0) is an unstable fixed 
point. The open chain (5J — J) on the other hand is a 
stable fixed point where the perturbation is described by 
the leading irrelevant operator d x cp(N)d x <j)(0) with scal- 
ing dimension of d = 2. Hence we expect a renormaliza- 
tion flow between the two fixed points as the temperature 
is lowered, and the temperature dependence of the impu- 
rity susceptibility as well as local properties will be de- 
scribed by a crossover function. Below a certain crossover 
temperature Tk this crossover function describes the be- 
havior of the stable fixed point (the open chain) while 
above Tk the system may exhibit a completely different 
behavior. The crossover temperature Tk itself is deter- 
mined by the initial coupling strength SJ 

lim T K -> 

SJ^O 

lim T K -> oo. (29) 

In other words, close to the unstable fixed point the 
crossover temperature is very small, indicating that we 
have to go to extremely low temperatures before we can 
expect to observe the behavior of the stable fixed point. 

A similar scenario holds for the impurity model with 
one external spin S t , 

H 2 = H + J'Sx-S/. (30) 

In this case, the periodic chain (J' = 0) is also the un- 
stable fixed point, while the open chain with, a decoupled 
singlet (J 7 — > oo) is the stable fixed pointO'EI 

As mentioned above, the renormalization flow has been 
confirmed for the energy corrections of individual eigen- 
states, but we now seek to extend this analysis to ther- 
modynamic properties, which will allow us to determine 
the crossover behavior and Tk directly. 

B. One weak link 

Our first task is to establish the renormalization behav- 
ior from a periodic chain fixed point to the open chain 



fixed point as a function of temperature. To examine 
the effective boundary condition on the spin-operators, 
it is instructive to look at the correlation functions at 
spin-sites close to the impurity. For periodic boundary 
conditions, the leading operator for the spin ^-operator 
is given by cos y/2nd> with scaling dimensions of d = 1/2 
according to Eq. (p7|). On the other hand, open boundary 
conditions restrict the allowed operators and the leading 
operator for Sf is found to be d x (f)(0) with scaling dimen- 
sion of d — 1 . Hence, the autocorrelation function at the 
impurity behaves differently, depending on the effective 
boundary condition 

wmwo» « { % ^°t b ' c ' < 31 > 

Therefore, a useful quantity to consider is the response 
Xi to a local magnetic field given by the Kubo formula 

Xi(T) = J 1/T (St(r)Sm) dr 

t^o J— logT periodic b.c. ,„„^ 

* \ const. + 0(T) open b.c. ^ ' 

In Fig. |^ we have presented the results of xi f° r different 
impurity strengths 5J on a logarithmic temperature scale. 
For the periodic chain (6J = 0) we clearly observe the 
logarithmic scaling, but for any finite 8J a turnover to a 
constant behavior is observed as T — ► 0, i.e. the behavior 
of the open chain. The turnover temperature Tk occurs 
at larger and larger values as we approach the stable fixed 
point [see Eq. (psj)]. An interesting aspect is that the 
curves actually cross: At very high temperatures the high 
temperature expansion always dictates a larger response 
for a weakened link, while at very low temperatures this 
relation is reversed by quantum mechanical effects and 
the renormalization flow. 




0.1 — 1 — ■ 1 — ■ 1 — ■ 1 — ■ 1 — 1 

-4 -3 -2 -1 

log(T) 

FIG. 2. The local response xi °f the spin Si 
next to a weak link for different coupling strengths 
SJ/ J = 0, 0.2, 0.4, 0.6, 0.8, 1. The crossover to the open chain 
behavior occurs at different temperatures Tk depending on 
the coupling strength. 
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We now turn to the true impurity susceptibility of 
Eq. (0) which is the experimentally more relevant quan- 
tity. Maybe the simplest non-trivial case to consider is 
the open chain SJ = J. At low temperatures the impurity 
susceptibility can be calculated from the leading irrele- 
vant local operator which is allowed in the Hamiltonian. 
This operator turns out to be (d x (t>(Q)) 2 which gives a 
constant impurity susceptibility with aJpgarithmic cor- 
rection similar to the pure susceptibilitytLZl as follows from 
a dimensional analysis. (In fact this operator can be ab- 
sorbed in the free Hamiltonian bw-a defining a velocity v 
that depends on the system sizella. An explicit calcula- 
tion of integrals over the correlation functions also comes 
to the same conclusion.) 



open 
^imp 



(T) 



const. + 0(1/ log(T/T )) 



(33) 



While it is possible to calculate the impurity susceptibil- 
ity numerically according to Eq. ( |2l| ) , the use of a second 
derivative in Eq. ( p2| ) causes large problems with the ac- 
curacy at lower temperatures since it involves taking the 
differences of large numbers. Luckily, the excess local 
susceptibility xiocai of the first site under a global mag- 
netic field turns out to-give a good estimate of the true 
impurity susceptibilityll3 



Xlocal 



d(sr, 

dB 



Xpurc oc Ximp + const. 



(34) 



where B is a global magnetic field and the constant is due 
to the alternating partliS. The results for this quantity are 
shown in Fig. |^ which are consistent with Eq. (p3|). The 
second derivative in Eq. (|2^) has a similar behavior in 
the intermediate temperature range, but is not accurate 
enough to extrapolate to the T — > limit as explained 
above. 
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FIG. 3. The impurity susceptibility for an open chain, 
which is approximated by the excess susceptibility at an 
open end Xiocai (dots). At moderately low temperatures our 
data- fit (solid line) is consistent with the behavior in Eq. (^). 
The fit is xiocai = 1.895 + 18.6/ log(T/T ), with log(T ) = 10. 



Now we are in the position to consider the impurity 
susceptibility of one weak link in the chain. By reducing 
SJ it is possible to tune the the system all the way from 
the open chain fixed point to the periodic chain. The 
operator (d x <f>(0)) 2 , which was responsible for the open 
chain impurity susceptibility is thereby reduced continu- 
ously. However, it is an entirely different operator corre- 
sponding to Sjy • Si which is responsible for the renormal- 
ization. This operator changes scaling dimension as we 
go from periodic boundary conditions towards the open 
chain 



Si oc 



sin V27T</> 
d x <l>(N)d x <K0) 



d= 1/2 
d = 2 



periodic b.c. 
open b.c. 



(35) 



Since we wish to study the effect of this renormalization, 
we choose to subtract the open contribution systemati- 
cally from Ximp, so we obtain exactly the part which will 
exhibit the crossover of the renormalization 



Xim P («5J) - SJ xZ7 = f(T/T K )/T K . 



(36) 



We see that this difference is zero at either fixed point. 
After subtracting SJx°m^ the impurity susceptibility 
comes only from the operator in Eq. (pa), which allows us 
to postulate the scaling form in Eq. Q36|). Below Tk the 
difference in Eq. ( p6| ) will asymptotically go to a constant 
for a very weak link SJ ~ J, which comes from the d = 2 
operator in Eq. ( |35| ) . Indeed we observe in Fig. ^ that this 
difference is negative and proportional to J — SJ ~ 1 /Tk 
and largely temperature independent as T — > 0. On the 
other hand as SJ becomes small enou gh, the behavior is 
much different: The expression in Eq. ( pq ) even decreases 
with SJ and the turn-over to constant behavior happens 
at much lower temperatures (outside the plot range). 




-0.04 



-3 



Log(T) 



FIG. 4. The part of the impurity susceptibility which is 
attributed to the effects of the weakened link. For a very weak 
link SJ ~ J the magnitude increases linearly with J — SJ and 
saturates quickly, while the magnitude decreases again as wc 
approach the periodic chain SJ <C J, but does not saturate as 
T -> 0. 
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Interestingly, this results in a highly nontrivial behav- 
ior as a function of T and Tk close to the unstable fixed 
point 



lim lim f(T/T K )/T K (xT K /T 2 oc SJ 2 -> 



lim lim f{T/T K )/T K cx 1/T K 



l/SJ 2 



(37) 



T/iis means that at T = a minute perturbation 6J re- 
sults in an extremely lavge negative Ximp 

l/SJ 2 (al- 
though this behavior occurs in an "unphysical" limit). 
The reason that the two limits do not commute is of 
course because one describes the behavior of the stable 
fixed point T <C Tk, while the other one describes the 
behavior at the unstable fixed point. While our numeri- 
cal results cannot show the entire crossover of Eq. (36), 
the increase below Tk for small J — SJ is clearly observed 
as well as the decrease and change of curvature above Tk 
for small SJ. Therefore, our data in Fig. [| supports the 
renormalization scenario and the nontrivial behavior of 



Eq. @. 



C. One external spin 



The model of an external spin antiferromagnetically 
coupled to the chain in Eq. (|l^) is maybe a little more 
exotic, i-ibAit, is still of great interest in a number of 
studiesacltJ. In Ref. [| it was first shown that the 
stable fixed point corresponds to open boundary con- 
ditions with, a decoupled singlet. This was confirmed 
numerically!], but more recently Dr. Liu postulated a 
completely different behavior using some non-local trans- 
formations on Fermion- fields which -mysteriously were re- 
arranged to form a solvable modeO. While we cannot 
trust or understand many of his calculations, the pre- 
dictions are in strong contrast to any previous expecta- 
tions and should be tested explicitly. In particular, he 
predicted the response of the impurity spin to a local 
magnetic field 



Xf(T) 



l/T 



(S* f (r)S* f (0))dT 



(38) 



to be proportional to T 5 / 2 at the Heisenberg point. We 
predict, however, that this response is described by the 
autocorrelation function of the leading operator for S f . 
By a symmetry analysis we find that for open boundary 
conditions this operator is given by d x 4>(0) with scaling 
dimension d = 1. The local response is therefore a con- 
stant as T — > with a linear term 



Xf (T) T -^2 const. 



0(T). 



(39) 



This also agrees with the findings in Ref. |^ which had 
similar reservations about Ref. yjj We now explicitly 
calculate Xf f° r several coupling strengths J' as shown 
in Fig. ||. Our data fits well to the predicted form and 



we can certainly rule out any T 5 / 2 behavior. Moreover, 
we find a scaling behavior which holds for all coupling 
strengths 



Xf (T) = g(T /Tk) /Tk- 



(40) 



A similar scaling relation was observed beforeu. Otar re- 
sults are consistent with previous numerical studiesQ, but 
distinctly larger in the low-temperature region. We at- 
tribute this to the finite size method used in Ref. ^, which 
becomes unreliable when the temperature falls below the 
finite size gap of the system. At large J' our results 
can be compared to that of two coupled spins forming 
a singlet. Note, that the response to a local magnetic 
field on one spin in a singlet is finite as T — > and pro- 
portional to 1/J' (and does not show activated behavior 
as a simple calculation shows). Therefore, our findings 
are completely consistent with the expectation that the 
impurity spin is locked into a singlet at the stable fixed 
point. 

log(T) 



1 - J'=1.5J 



0.10 



0.05 - 



0.00 




FIG. 5. The upper 
of the impurity spin Xf 
J' / J = 0,0.3,0.6,0.9,1.2. 



1 
log(T/T K ) 

figure shows the local response 
for different coupling strengths 
1.5. The lower figure shows the 



data collapse according to Eq. (^). In the inset the corre- 
sponding crossover temperature Tk is shown, which is ex- 
pected to be T K oc J' 2 for J' < J and T K oc J' for J' > J. 



IV. DISCUSSION 

To accurately calculate impurity properties we have to 
determine not only the largest eigenvalue of the transfer 
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matrix, Tm, but also the corresponding eigenvectors to 
high accuracy. To estimate the error of the impurity 
properties is difficult. Errors come both from the finite 
Trotter number, M, and from the finite number of states, 
m, in the DMRG. The scaling of pure properties with M 
and to usually turn out to be simpler than the scaling of 
impurity properties, which show a less clear form of the 
errors. 

We have used the value (3/M = 0.05 for all calculations 
presented in this article. With this value, the error due 
to the finite M should be small, which is also confirmed 
by test-runs. 

We have tested the error due to the finite Trotter num- 
ber M by doing separate DMRG runs for different values 
of (3/M . For the pure case at moderate temperatures we 
find that the eigenvalue A scale as 1/M 2 , as is expected, 
while at the lowest temperatures, the error due to the 
finite m is larger making it difficult to see the expected 
1/M 2 scaling. For the impure case, the convergence of 
Aimp with M is more complicated, but the overall scaling 
is however still roughly \/M 2 . 

We have also tested the convergence with the number 
of basis states m. For both the pure and the impure cases 
we find a rapid convergence with increasing m. We have 
used a maximum of m = 65 for the calculations on the 
weakened link impurity and m = 38 in the external spin 
case. We found however no noticeable difference between 
777 = 38 and m = 65 down to T = 0.02 in a test run for 
the external spin. 

The truncation error, i.e. 1 — ^2 i=1 Wi, where Wi are 
the largest eigenvalues of the density matrix, is less than 
about 10~ 5 for m = 65 at the lowest temperatures. Note 
that the truncation error is determined during the pure 
sweep in which also the projection operators and target 
states are determined. It could thus be used as an esti- 
mate of an upper limit of the error of the pure properties, 
but it is difficult to say how good estimate it yields for 
the impurity properties. 

To test our results we have also done Quantum Monte 
Carlo (QMC) simulations for a few temperatures and 
couplings. The DMRG data was well within the error 
bars of the QMC results. 

Local properties converge much faster with m than 
the impurity susceptibility. This fact might be explained 
by the difficulty to numerically take a second derivative, 
since we have to subtract two large numbers to find SF 
in Eq. (^2|). This is however not the case for local prop- 
erties, since we know that the local magnetization is zero 
in the absence of a magnetic field. Another source of 
inaccuracy in Ai mp comes from the error of the target 
states. Let us assume that the target state is determined 
up to some error e; IV') exact = IV^dmrg + k). Since 
the target state is, to numerical accuracy, an eigenstate 
of Tm, the eigenvalue will be determined to order e 2 . 
Expectation values of other operators, for example the 
impurity transfer matrix, will however in general only 
be accurate to order e. The local properties don't seem 
to suffer too much from this effect, the reason might be 



that there is some cancellation of errors in the quotient 

While the accuracy of the second derivative is good 
enough for the pure susceptibility down to about T — 
0.01, we cannot trust the impurity susceptibility below 
temperatures roughly an order of magnitude larger. Lo- 
cal impurity properties on the other hand seem to be well 
represented down to about T = 0.02. 

In summary we have shown that the transfer matrix 
DMRG is a useful method for calculating finite temper- 
ature impurity properties of a spin chain in the thermo- 
dynamic limit. We have considered two impurity mod- 
els: One weakened link and one external spin, but the 
method can be applied to other impurity configurations 
and electron systems. We find that the local response of 
the spin next to a weakened link always crosses over to a 
constant below some Tk, i.e. to the behavior of the open 
chain fixed point. According to our calculations, the im- 
purity susceptibility shows an exotic cross-over behavior 
with non-commuting limits. For the external spin impu- 
rity, we have found that the data for the local response 
shows the expected cross-over to open chain behavior as 
the temperature is lowered. The response has the scal- 
ing form in Eq. ( |40| ) and we can explicitly show the data 
collapse and determine Tk (Fig. ||). 
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